# forecasts.R
# This file is the batch file to run the BVAR analysis and conditional
# forecasting for the Levant model
#
# Patrick T. Brandt
# pbrandt@utdallas.edu
#
# September 2004
# Revised: July 2005 to work with latest version of MSBVAR
#

# Load libraries 
library(foreign)  # Functionts to read the Stata file of the raw data
library(GenKern)  # Library that can draw 3d and 2d correlated
                  # k.p.d.f.s for the mountain plots

library(MSBVAR)   # Brandt's package for Bayesian VAR models

# Read the data and set up as a time series
data <- read.dta("levant.weekly.79-03.dta") 
attach(data)

# Set up KEDS data
KEDS.data <- ts(cbind(a2i,a2p,i2a,p2a,i2p,p2i),
                start=c(1979,15),
                freq=52,
                names=c("A2I","A2P","I2A","P2A","I2P","P2I"))

# Select the sample we want to use.
KEDS <- window(KEDS.data, end=c(1988,50))


################################################################################
# Estimate the BVAR models 
################################################################################

# Fit a flat prior model
KEDS.BVAR.flat <- szbvar(KEDS, p=6, z=NULL, lambda0=1,
                         lambda1=1, lambda3=1, lambda4=1, lambda5=0,
                         mu5=0, mu6=0, nu=0, qm=4, prior=2,
                         posterior.fit=F)

# Reference prior model -- Normal-IW prior pdf
KEDS.BVAR.informed <- szbvar(KEDS, p=6, z=NULL, lambda0=0.6,
                             lambda1=0.1, lambda3=2, lambda4=0.5, lambda5=0,
                             mu5=0, mu6=0, nu=ncol(KEDS)+1, qm=4, prior=0,
                             posterior.fit=F)

# Set up conditional forecast matrix conditions
nsteps <- 12
a2i.condition <- rep(mean(KEDS[,1]) + sqrt(var(KEDS[,1])) , nsteps)

yhat<-matrix(c(a2i.condition,rep(0, nsteps*5)), ncol=6)

# Set the random number seed so we can replicate the results.
set.seed(11023)

# Conditional forecasts
conditional.forcs.ref <- hc.forecast.var(KEDS.BVAR.informed, yhat, nsteps,
                            burnin=3000, gibbs=5000, exog=NULL)

conditional.forcs.flat <- hc.forecast.var(KEDS.BVAR.flat, yhat, nsteps,
                             burnin=3000, gibbs=5000, exog=NULL)

# Unconditional forecasts
unconditional.forcs.ref <-uc.forecast.var(KEDS.BVAR.informed, nsteps,
                                          burnin=3000, gibbs=5000)

unconditional.forcs.flat <- uc.forecast.var(KEDS.BVAR.flat, nsteps,
                                            burnin=3000, gibbs=5000)


save.image("forecasts.RData", compress=T)

# Set-up and plot the unconditional and conditional forecasts.  This
# code pulls for the forecasts for I2P and P2I and puts them into the
# appropriate array for the figures we want to generate.
uc.flat <- NULL
hc.flat <- NULL
uc.ref <- NULL
hc.ref <- NULL

uc.flat$forecast <- unconditional.forcs.flat$forecast[,,5:6]
hc.flat$forecast <- conditional.forcs.flat$forecast[,,5:6]
uc.ref$forecast <- unconditional.forcs.ref$forecast[,,5:6]
hc.ref$forecast <- conditional.forcs.ref$forecast[,,5:6]

actual <- window(KEDS.data, start=c(1988, 51), end=c(1989,10))[,5:6]

# Build a special function for the forecast plots so that we can
# present the two forecasts and the actual data.  This is a
# modification of the plot.var.forecasts() function in the MSBVAR
# code.

# Here the fcasts1 and fcasts2 are from the Gibbs sampled forecasts.
# The fcasts3 are the actual levels

plot.var.forecasts3 <- function(fcasts1, fcasts2=NULL, actual=NULL,
                                varnames=NULL,
                                start=c(0,1), 
                                freq=1, probs=c(0.05,0.95),
                                compare.level=NULL, ylab=NULL, ...)
  
  { # compute quantities for ecdf of forecast matrix 1
    m <- dim(fcasts1$forecast)[3]
    h <- dim(fcasts1$forecast)[2]
    iters <- dim(fcasts1$forecast)[1]
    fcast1.summary <- array(apply(fcasts1$forecast, 3, forc.ecdf, probs=probs), c(h,3,m))

    # Now do the same for forecast 2 if non-NULL
    if (is.null(fcasts2)==FALSE)
      { fcast2.summary <- array(apply(fcasts2$forecast, 3,
                                      forc.ecdf, probs=probs),
                                c(h,3,m))
      }
    
    par(las=1, mar=c(1,2,2.5,1))
    for(i in 1:m)
      { forc1.ci <- ts(fcast1.summary[,,i], start=start)
        
        if(is.null(fcasts2)==FALSE)
          { forc2.ci <- ts(fcast2.summary[,,i], start=start) }

        if(is.null(actual)==FALSE)
          { forc3.ci <- ts(actual[,i], start=start) }
        
        ylim <- c(floor(min(c(forc1.ci,forc2.ci,forc3.ci))),
                  ceiling(max(c(forc1.ci,forc2.ci,forc3.ci,compare.level))))
        
        ts.plot(forc1.ci, forc2.ci, forc3.ci,
                gpars=list(lty=c(1,1,1,2,2,2,4),
                   col=c(rep("black",3),rep("black",3), "black"),
                   ylim=ylim, xlab="",axes=FALSE, ... ))
        axis(2,c(floor(min(c(forc1.ci,forc2.ci))),
                 ceiling(max(c(forc1.ci,forc2.ci)))))
        mtext(varnames[i],side=3,line=1)

        box();
        if(i==1) { mtext(ylab, side=2, line=3, at=c(1.5*mean(ylim))) }
        abline(h=0)

        # put in the comparison level if one is provided
        if (is.null(compare.level)==FALSE)
            { abline(h=compare.level[i], lty=c(2)) }
       
      }
#    par(oldpar)
  }


# Conditional and unconditional forecasts for I2P and P2I for the BVAR
# paper.  These are figures 3 and 4 in the paper.

pdf(file="cndforecasts.pdf", width=7.5, height=5) 
par(mfrow=c(2,2), omi=c(0.25,0.5,0.25,0.25)) 
plot.var.forecasts(uc.flat,hc.flat,
                    probs=c(0.16, 0.84),
                    varnames=c("I2P", "P2I"),
                    compare.level=KEDS[nrow(KEDS),5:6], lwd=2)
plot.var.forecasts(hc.ref,hc.flat,
                   probs=c(0.16, 0.84),
                   varnames=c("I2P", "P2I"),
                   compare.level=KEDS[nrow(KEDS),5:6], lwd=2)
dev.off()

pdf(file="uncndforecasts.pdf", width=6, height=5) 
par(mfrow=c(1,2), omi=c(0.25,0.5,0.25,0.25)) 
plot.var.forecasts3(uc.flat,uc.ref,actual,
                   probs=c(0.16, 0.84),
                   varnames=c("I2P", "P2I"),
                   compare.level=KEDS[nrow(KEDS),5:6], lwd=2)
dev.off()



# Now make the mountain plot for the forecasts in week 12.  This is
# the final figure in the BVAR / PA paper.


# Do the kernel smoothing for the reference prior and flat prior
# forecasts.

ref <- KernSur(conditional.forcs.ref$forecast[,12,5],
               conditional.forcs.ref$forecast[,12,6],
               xgridsize=40, ygridsize=40)
flat <- KernSur(conditional.forcs.flat$forecast[,12,5],
                conditional.forcs.flat$forecast[,12,6],
                xgridsize=40, ygridsize=40)

ref.ip <- KernSec(conditional.forcs.ref$forecast[,12,5])
ref.pi <- KernSec(conditional.forcs.ref$forecast[,12,6])

flat.ip <- KernSec(conditional.forcs.flat$forecast[,12,5])
flat.pi <- KernSec(conditional.forcs.flat$forecast[,12,6])


pdf(file="mountain.pdf", width=6, height=6)
par(mfcol=c(2,2))
par(mai=c(0.75,0.75,0.2,0.2))

plot(ref.ip$xords, ref.ip$yden,
     ylim=range(ref.ip$yden, flat.ip$yden),
     xlim=range(ref.ip$xords, flat.ip$xords),
     type="l", xlab="I2P", ylab="Density", lwd=2)
lines(flat.ip$xords, flat.ip$yden, lty=2, lwd=2)

plot(ref.pi$xords, ref.pi$yden,
     ylim=range(ref.pi$yden, flat.pi$yden),
     xlim=range(ref.pi$xords, flat.pi$xords),
     type="l", xlab="P2I", ylab="Density", lwd=2)
lines(flat.pi$xords, flat.pi$yden, lty=2, lwd=2)

contour(ref$yords, ref$xords, ref$zden, lty=1, lwd=1,
   xlim=range(-60,50), ylim=range(-60,50),
   zlim=range(ref$zden,flat$zden), xlab="", ylab="", drawlabels=F)
abline(h=0)
abline(v=0)
par(new=T)

contour(flat$yords, flat$xords, flat$zden, lty=c(2), lwd=1,
        xlim=range(-60,50), ylim=range(-60,50),
        zlim=range(ref$zden,flat$zden), xlab="P2I", ylab="I2P", drawlabels=F)

par(mai=c(0.5,0.25,0.1,0.25), cex=0.4)

fill <- matrix("grey", nr = nrow(flat$zden) - 1, nc = ncol(flat$zden)
               - 1)
fill[, i2 <- c(1, ncol(fill))] <- "gray"
fill[i1 <- c(1, nrow(fill)), ] <- "gray"

persp(flat$xords, flat$yords, flat$zden,
      col=fill, shade=0.5,
      border="black",
      xlim=range(ref$xords,flat$xords),
      ylim=range(ref$yords,flat$yords),
      zlim=range(ref$zden,flat$zden),
      xlab="I2P", ylab="P2I", zlab="",
      ticktype="detailed", theta=225)
par(new=T)

persp(ref$xords, ref$yords, ref$zden,
      col="transparent",
      border="black",
      xlim=range(ref$xords,flat$xords),
      ylim=range(ref$yords, flat$yords),
      zlim=range(ref$zden,flat$zden),
      xlab="", ylab="", zlab="",
      shade=0.6, ticktype="detailed",  theta=225)


dev.off()


# End of file








